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We propose a simple model that elucidates the generation of halo velocity bias. The fluid equation 
approximation is often adopted in modelling the evolution of the halo density field. In this approach, 
halos are often taken to be point particles even though in reality they are finite-sized objects. In 
this paper, we generalize the fluid equation approximation to halos to include the finite extent of 
halos by taking into account the halo profile. We compute the perturbation of the halo density and 
velocity field to second order and find that the profile correction gives rise to correction terms 
in Fourier space. These corrections are more important for velocity than for density. In particular, 
the profile correction generates correction term in the velocity bias and the correction terms do 
not decay away in the long term limit, but it is not constant. We model the halo profile evolution 
using the spherical collapse model. We also measure the evolution of proto-halo profile at various 
redshifts numerically. We find that the spherical collapse model gives a reasonable description of the 
numerical profile evolution. Static halo profile is often adopted in modelling halos in theories such 
as the excursion set theory. Our work highlights the importance of including the profile evolution 
in the calculations. 


I. INTRODUCTION 

Peculiar velocity can be a powerful probe of cosmol¬ 
ogy. On one hand, peculiar velocity causes redshift space 
distortion [T], thus to get cosmological information from 
galaxy surveys, one needs to model the peculiar veloc¬ 
ity as well. On the other hand, redshift space distor¬ 
tion causes anisotropy, which gives rise to higher order 
multipoles in the correlation function/power spectrum 
[2], and hence a useful signal. Peculiar velocity can also 
help distinguish general relativity from modified gravity, 
e.g. Hi]. The redshift space distortion has been used to 
constrain the growth rate and testing gravity from some 
recent galaxy surveys [Mil]. In future the peculiar ve¬ 
locity surveys can also be fruitful [12) . 

As in galaxy surveys, only galaxies are observable, not 
the dark matter, it is important to understand if the ve¬ 
locity of the galaxy is biased with respect to that of the 
underlying dark matter or not. As galaxies are hosted in 
halos, and halos are simpler than galaxies because they 
are only governed by the gravitational physics, in this 
paper, we take studying the halo bias as a step towards 
understanding the galaxy bias. Recently there have been 
some indications from the measurements using halos at 
low redshifts that velocity bias may be non-negligible, 
although the quantitative measurement is still hard [T^ 
[n]. Velocity measurement is nontrivial because it re¬ 
quires the velocity field of the tracers to be weighted 
by volume, and it is easy to mistakenly get the density 
weighted velocity, i.e. momentum, instead of velocity 
m- When the number density of the tracers is low, it 
suffers from numerical sampling artifact, see e.g. [13 m- 
121) . One way out is to use momentum instead, e.g. in 
[mEg. Unlike the velocity field, however, there is an 
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additional complication that in momentum the galaxy 
density bias is involved as well. 

On the theory side, in the usual fluid approximation 
for dark matter and galaxy, even if the initial velocity 
field of the galaxy differs from that of the dark matter, 
i.e. there is initial velocity bias, large scale gravitational 
evolution will naturally drive the galaxy velocity field to 
that of the dark matter [IIE3]- On the other hand, 
the peak model predicts that the velocity bias persists 
and remains constant at late time [21] . This result seems 
to be favoured by the recent simulation results, which 
suggest that the halo velocity bias at late time is non- 
negligible at fc ^ O.I5Mpc“^/i. Ref. m argued that 
the force on the halo has to be “biased” in order for the 
coupled-fluid approach to agree with peak model result, 
although no further justification was given. Ref. [21] tried 
to derive the peak theory results using the distribution 
function approach. Here we take a different approach. 
A halo is a composite object consisting of a collection 
of particles. The position of the halo is defined by the 
position of the center of mass (CM) of its constituent 
particles. Thus the force acting on the CM position of the 
halo should be averaged over its constituent particles. In 
this way, we give a physical origin for the “biased” force 
on the halo. We will show that the halo profile correction 
naturally gives rise to the leading correction to the 
velocity bias and it does not decay away. 

We note that our approach also has rather different in¬ 
terpretation for the generation of velocity bias from that 
in In peak theory, although the smoothing 

window is an important ingredient, the window function 
is usually assumed to be static. Sometimes, the attention 
is focused on the discrete peak “points”, which have the 
same velocity as the dark matter locally, it was argued 
that the velocity bias is a “statistical” effect. In our 
model, the velocity bias physically arises from the fact 
that halos are finite-sized objects, not point particles. 
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and it also highlights the dynamical nature of window. 

On the other hand, our approach may not be mutually 
exclusive with the peak model approach. In the mod¬ 
elling of halos starting from the Lagrangian space, one 
defines halos with window function and the smoothing 
scale is fixed to be the Lagrangian szie when they are 
transformed to the Eulerian space. Even if our velocity 
bias contribution is not the dominant one seen in simu¬ 
lations, the profile correction effects should be taken into 
account in the calculations as well. 

This paper is organized as follows. As we will show 
that the halo prohle gives rise to the velocity bias cor¬ 
rection, to set the stage, we will first review the evolu¬ 
tion of the halo profile using the spherical collapse model 
in Sec. |II A| and the numerical halo profile is measured 
from simulation and compared with the spherical col¬ 
lapse model in Sec. \\rB\ In Sec. [IITAI we compute the 
correction to the linear velocity and density bias due to 
the halo profile, and the second order corrections are pre¬ 
sented in Sec. IIIIBI We conclude in Sec. ||Vl 


II. HALO PROFILE EVOLUTION 

Halo prohle is often used in the context of halo model 
for modelling the dark matter power spectrum [271130] . 
In this case, the virialized halo prohle, such as the NEW 
prohle [3T] is often used. However, we will follow the 
proto-halo from its infancy to the hnal virialized stage in 
modelling the bias evolution. To this end, we will hrst 
review the evolution of a halo using the spherical col¬ 
lapse (SC) model. We will then construct proto-halos at 
various redshifts and measure the prohle evolution in nu¬ 
merical simulations. The results are compared with the 
SC model. To our knowledge, this is the hrst systematic 
measurements of the proto-halo prohle evolution. 


A. Profile evolution from SC model 


A simple analytic model for halo evolution is given 
by the SC model [32] (see also [33H35] '). Suppose that 
the initial huctuations are spherically symmetric about 
some point in position space. To avoid shell crossing, 
we assume that the radial prohle is non-increasing as 
the distance from the center increases. We will consider 
the matter dominated universe as the resultant equation 
can be integrated analytically, and also the more realistic 
ACDM model. Under the Newtonian approximation, the 
equation of motion for a mass shell at a distance r from 
the center is given by 

d^r GM(r) A 

where both r and t are the physical distance and time, 
and G is the gravitational constant, M{r) is the total 
mass inside the mass shell and A is the cosmological con¬ 
stant. 


Integrating Eq. [^ once, we obtain the hrst integral of 
motion 


1 /dr\'^ GM Ar^ 
2\dt) r 6 


( 2 ) 


where the total energy U is a constant of integration. 

We will solve Eq. [^ numerically when A 0. When 
A = 0, Eq. [^can be further integrated analytically, and 
the solution can be expressed in the form of a cycloid 
solution 


r = A(l — COS0), (3) 

t + T = B{0 -sin9), (4) 

= GMB^, (5) 

where A, B and T are constants. The parameter 6, also 

called the development angle, runs from 0 to 27r. When 9 
is close to 0, the overdensity inside the mass shell is small, 
the mass shell essentially follows the Hubble expansion. 
The mass shell reaches maximum at 9 = n. Beyond 
that the mass shell overcomes the Hubble expansion and 
turns around. At 9 = 2 tt, the shell collapses to a point 
according to Eq. [^ However, it was argued that during 
the rapid infall the potential varies quickly, the particles 
no longer follow the energy conserving orbits, instead the 
energy available to the particles is widened, and the sys¬ 
tem reaches virial equilibrium [3S1[37]. This procoess is 
called “violent relaxation” [36] . From virial equilibrium, 
one finds that the virial size is related to as 



The virial size is often given in terms of the virial den¬ 
sity, Av, as 



where M is the mass of the halo and pm is the comoving 
density of matter. For EdS universe, Av is equal to 178. 
When A 7 ^ 0, a fitting formula for the virial density, Av, is 
given in [SS]. For the flat ACDM with Dm = 0.25 adopted 
in this paper, Av for a halo virailzed at z = 0 is 380. In 
practice other values of Av are often adopted, such as 200 
and 500. We will use Av = 500 as we will see later on 
it gives a good description of our simulation data. Note 
that for non-EdS universe, in Eq. [^ the critical density 
is often used instead of Pm to define r^. 

To solve Eq. [^ the initial conditions that the ini¬ 
tial overdensity is obtained by extrapolating the collapse 
threshold from the present time to the initial time using 
the linear growth factor and zero initial peculiar velocity 
is often assumed. However, the zero initial peculiar ve¬ 
locity condition excites both the growing mode and the 
decaying mode. It can be shown that the linear ampli¬ 
tude of the perturbation is reduced by a factor of 3/5 
[321 [Ml 12S] • This is equivalent to setting up the initial 
condition incorrectly and transient effects are induced. 
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To get the right final amplitude, one quick fix is to in¬ 
crease the initial perturbation by a factor of 5/3 to com¬ 
pensate the loss to the decaying mode. As we set up the 
initial conditions at not very high redshifts, the transient 
effects are not negligible. A better approach is to set 
the initial peculiar velocity such that the decaying mode 
vanishes. Thus we will use the initial conditions 
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(8) 

r 

VdTTpin^ 

/ 

(1 + 5*)5’ 

1 - \ 

r* = 1 

-T-)’ 

(9) 


In this paper, we use to denote a quantity at some 
initial time. Thus a*, i/* and (5, are the scale factor, 
Hubble parameter and the average density contrast inside 
the spherical shell at the initial time. 

Using these initial conditions, we can write the coeSi- 
cents A and H in Eq. and as 


A = 




^m(l + 


2 11 ^( 1 + 


2 ’ 




^m(l + ^*) 


ij* + s,) - (1 - iiy 


( 10 ) 

(11) 


where is the density parameter of matter at the initial 
time. We first note that A and hence r is proportional 
to r*. On the other hand, the collapse history given by t 
is independent of r*, and it depends only on the matter 
inside through UJ/ and (5*. 

In Eq. after dividing by r», one can easily see that 
the collapse history is independent of r* in ACDM model. 
Thus given and ft a, the collapse history depends only 
on 6^. 

In Eig.[^ we show the evolution of the profile as a func¬ 
tion of a for three different cosmological models: EdS, 
Open CDM with Om = 0.25 and ACDM with Om = 0.25 
and Da = 0.75. Note that in this paper ACDM always 
refers to this model. The mass of the halo is chosen to 
be 2 X 10^^ The collapse threshold at z = 0 is 

set to be i5c = 1.68 and extrapolated to the initial time 
using the linear growth factor for the corresponding cos¬ 
mology. The initial conditions are set using Eq. and 
1^ Nonetheless, we still find that we need to choose a* 
to be sufficiently small (a* = 0.01 here) to reduce the 
effects of transients. For example when a* = 0.02 is 
chosen instead, we find that the collapse epochs are in¬ 
creased by a few per cent compared to the ones shown. 
We emphasize that this is because we set up the initial 
conditions using linear theory, and the transients can be 
further suppressed using higher order perturbation the¬ 
ory. This is analogous to setting up initial conditions 
in simulations using 2LPT [sasoi- Although these cos¬ 
mologies are rather different, the final collapse epochs are 
very similar as long as the correct linear growth factor is 
used to set 5*. In other words the collapse threshold is 
insensitive to the cosmological model [H]. However, the 
intermediate stages of collapse are quite different among 


these models. Thus this suggests that we need to solve 
the model explicitly in order to follow the evolntion of 
the halo profile accurately. We also indicate in Fig. 
the virial size after the collapse. For EdS and OCDM, it 
is computed using Eq. while for ACDM we use Eq. 
with Av = 500. 



FIG. 1. The evolution of the halo profile in three different 
cosmological models EdS (blue), OCDM with Dm = 0.25 (red) 
and ACDM with Um = 0.25 and Um = 0.75 (green). The 
virial size is indicated as a dot after the collapse. 

To illustrate the evolution of the profile, we further 
assume that the initial halo profile is described by a top- 
hat profile, and in Fourier space, it reads 

3 

IUth(k) = ^ (sins; — K cos k). (12) 

Top-hat profile is a good approximation at high redshifts. 
In SC with top-hat perturbation, the top-hat shape is 
preserved during evolution. The only part that changes 
is the width of the window. On the other hand, the 
Eulerian virialized spherical halo profile is well described 
by the NEW profile |3I]. This means that a top-hat 
window of perturbation cannot evolve to the NFW-like 
profile. We will see in Sec. |HB| that the halo profile 
measured from simulation goes from one resembling top- 
hat to an NFW-like profile as the redshift decreases. 

In Fig.[^ we show the evolution of W^Yiikx) for a series 
of values of the comoving size x at different time a. We 
start with x* = 4.07Mpc/i“^ at z* = 99, and the size of 
the spherical shell is then evolved according to Eq. We 
have adopted a flat ACDM model with Da = 0.75, and 
5c = 1.68. Note that although the physical size r first 
expands and then collapses as in Fig. the comoving 
size X is always decreasing as the effect of expansion is 
removed. At a = 1, the shell has not fully collpased yet. 
As we will see later on, we are mainly interested in the low 
k part of the window, thus to a very good approximation, 
the window is essentially 1 up to fc ^ 3Mpc“^ h at the 
present time. In this plot we have not substituted the 
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halo size at o = 1 with the virial size. Most of time in the 
paper, the sudden change during virialization does not 
matter as it occurs almost instantaneous for our purpose. 
From now on, the SC model refers the one obtained by 
evolving an initially top-hat perturbation using Eq. 



FIG. 2. The evolution of the top-hat window, WTiiikx) for 
a suite of scale factor a. The comoving size x of the spherical 
shell is computed using the SC model. 


B. Measurement of the profile evolution from 
simulations 

In this section, we will measure the evolution of halo 
profile from A^-body simulation. Starting from Eulerian 
halos, such as those at z = 0, we trace the particles in 
the Eulerian halo back in time to construct the proto¬ 
halos at earlier times. The position of the proto-halo at 
redshift z is defined by the CM of its constituent particles 
at redshift z. We will consider proto-halos at various 
redshifts. 

Before presenting the numerical results we would like 
to first outline the details of the iV-body simulation used 
here. We shall use the Oriana and Carmen simula¬ 
tions in the LasDamas project. In these simulations, 
a flat ACDM model with the cosmological parameters, 
rini = 0.25, Ha = 0.75 and as — 0.8 are adopted. The 
transfer function is output from CMBEAST |32]- The 
initial conditions are Gaussian with spectral index be¬ 
ing 1. The initial displacement fields are created using 
2LPT [ID] at z = 49. The simulations are evolved using 
the code Gadget2 [33 . In the Oriana simulations, there 
are 1280^ particles in a cubic box of size 2400 Mpc/i“^, 
while for Carmen simulation there are 1120^ particles in 
a box of size lOOOMpc/i”^. Thus the particle masses 
are 4.57 x 10^^ and 4.94 x 10^° for Oriana and 

Carmen respectively. We shall use five realizations for 
Oriana and seven for Carmen. The halos are obtained 
using Friend-of-Friend halo finder. For Oriana, the link¬ 


ing length b = 0.156 is used, while b = 0.2 for Carmen. 
To resolve the halo better, we use halos with at least 150 
particles. Although the Carmen simulations have better 
mass resolution than Oriana, we find that their results 
are quite similar. To avoid redundancy, most of the time, 
we only show results from Oriana. 

In Fig. we show the halo profile at redshifts, z = 0, 
0.34, 0.73, 0.97, 1.5, and 49. The halo profile is obtained 
by stacking the halos in the same mass bin together and 
spherically averaged to get the spherically symmetric pro¬ 
file. In this plot, the Eulerian halo is at z = 0 and the 
proto-halos at higher redshifts are constructed from the 
Eulerian ones. We find that when the size of halos of 
different masses normalized by their corresponding virial 
size, they coincide well with each other. We computed 
Tv using Av = 500 although this is immaterial to our 
purpose here. 

Note that for redshift z = 0, the halos are in fact in 
the Eulerian space. The virialized spherical Eulerian halo 
profile at low redshift is well fitted by the NFW pro¬ 
file. However, the halos used to construct this profile are 
carefully selected, see e.g. (DU ED- These halos are con¬ 
structed using spherical overdensity finder and they are 
chosen to be spherically symmetric and in a relaxed state 
without signatures of recent mergers. Here we use all 
the halos obtained from the halo finder without further 
screening. We find that our Eulerian profile is reasonably 
well fitted by the NFW profile, but we also find that the 
profile close to the virial radius drops faster than r“^, the 
scaling of the NFW profile near virial radius. Using the 
NFW profile, one can show that the Eulerian profile is 
approximately universal in the variable r/xy for different 
masses because of the fact that the concentration only 
weakly depends on the mass of the halo m- 

The deviation of the proto-halo profile from the NFW 
profile increases as the redshift increases. At z = 49, 
the proto-halo profile corresponds to the one in the ini¬ 
tial condition of the simulation. Theoretically, the La- 
grangian profile is often assumed to be a top-hat. In 
[45] ■ it is found to be in between a Gasussian and a top- 
hat. More precisely, in Fourier space, the Lagrangian 
profile is well fitted by a product of a Gaussian and a 
top-hat window. We note that as z increases, there are 
large deviations in the profile at small r among different 
halo masses. 

As there is no universal halo profile that can fit the 
proto-halo profile at various redshifts well, we shall use 
the numerical profile directly. The profile is Fourier 
transformed numerically. In the case of the NFW profile, 
the integration is cut-off at the virial radius Vy |30j . In 
our case, for proto-halos at intermediate redshifts, i.e. in 
between z = 0 and z = 49, it is not clear what the cut-off 
size should be. Nonetheless, our profile drops rapidly for 
r greater than a few Xy, we can take r to be infinity, and 
the results are unaffected. As the mass of the proto-halo 
is conserved, we expect the low-fc part of the profile at 
different redshifts to be the same. Numerically, however, 
this is not always achieved. In each set of simulation, we 
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FIG. 3. The halo profile at j* = 0, 0.34, 0.73, 0.97, 1.5, and 49. The Eulerian profile is at z = 0, and the proto-halo profiles 
are obtained by tracing the particles in the Eulerian halo to higher redshifts, Zl. Results from halos of mass ranging from 
6.9 X to 1.3 X 10^® are shown. The density is normalized with respect to the mean comoving density of 

matter, pm and the radial distance r is normalized with respect to the Eulerian virial size of the halo, r^. The results are from 
the Oriana simulations. 


find that the fractional deviation of the p{k = 0) across 
redshifts decreases as the mass of the halo increases. For 
example, for the lowest mass halo used for Oriana, the 
fractional deviation of p{k = 0) is within 10%. This is 
one of the indications that we should use halos with large 
number of particles. From now on, we simply normalize 
the profile so that it is 1 at low k. In Fig. we show 
the Fourier transform of the halo profile for a selection 
of halo masses. As the redshift zl increases, the size of 
the proto-halo increases, the window in Fourier space de¬ 
creases and the low-A: plateau shrinks. We also note that 
there are oscillations in the Fourier transform of the win¬ 
dow. It is more visible as the mass of the halo increases 
because the oscillations are pushed to lower k. They are 
also more prominant as zl increases because the halo 
profile is more top-hat-like, and hence the wiggles are 
stronger. 


To show the scale-dependence of the window function, 
in Fig. we plot |kF — 1| obtained from simulations and 
the SC model. Absolute value is taken because the k'^ cor¬ 
rection is negative. We have introduced the time variable 
y = In D, where D is the linear growth factor (defined by 
Eq. 16). It is normalized such that ?/ = 0 at z* = 49. As 
we will see in next section, this time variable is conve¬ 
nient. First at y = 0, SC model agrees with data well ex¬ 


cept for the highest mass bin shown (3.7 x 10^® 
Unfortunately, we have no simulation data available in 
the range 0 < y < 3, although we expect that the over¬ 
density is still in the expansion stage (y < 3), the SC 
model should work reasonably well. When the region 
turns around and collapses, we expect the SC model to 
fail to describe the simulation data accurately. In fact, 
during the turn-around and collapse phase, the SC results 
are larger than the simulation data. We also note that for 
various values of k shown, the agreement between the SC 
model and the data is qualitatively smiliar. In Fig. we 
show the corresponding results obtained using the Car¬ 
men simulations. In this plot, the Eulerian halos are at 
z = 0 and proto-halos are constructed at z = 0.13, 0.52, 
0.97 and 49. Although Carmen has better mass resolu¬ 
tion, the results are quite similar to those obtained from 
Oriana. Overall, the agreement between the simulation 
results and SC model is reasonable. 

The SC model only works qualitatively at late time. 
After all, in SC the halo profile shape does not change, 
but for real halos the halo profile shape does change. In 
[46] , a modified SC model that tries to overcome the jump 
at the final virialization stage was proposed. The model 
joins smoothly to the final virial scale at the expense of 
two additional free parameters. The modified model is 
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2^=0, 2l=0 2j;=0, 2l=0.34 2j;=0, Zl =0.73 





FIG. 4. Same as Fig. except in Fourier space. The halo profile is normalized such that it approaches 1 as A: tends to 0. 


M=6.9e + 13 


M=2.0e + 14 


M = 6.1e + 14 


M=3.7e + 15 




FIG. 5. The function |IF — 1| as a function of y, obtained from simulations (filled circles) and SG model (solid line). Results 
from four mass bins of mass 6.9 x 10^^, 2.0 x 6.1 x 10^“^, and 3.7 x 10^® are shown (from left to right). For each 

mass bin, \W{k,y) — 1| at six different k’s are plotted. The data is from Oriana. 


valid only when the density is high as it is an expansion in 1 /6. Nonetheless, using these additional parameters, one 
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M=l.le + 13 


M=4.6e+13 


M=2.0e+14 


M=5.8e+14 







FIG. 6. Similar to Fig. except for Carmen. 


may get a profile evolution history, especially the part 
from turn-around to collapse, to agree with the simula¬ 
tion results better. Another potential way to improve the 
modelling is to use ellipsoidal model [331 ISSl EZl EE] • For 
example, the halo mass function motivated by the ellip¬ 
soidal collapse improves the agreement with simulation 
[33] compared to the spherical Press-Schechter one. The 
halo collapse threshold is also better modelled by the el¬ 
lipsoidal collapse model [50]. However, as the halo profile 
considered here is spherically averaged, one still need to 
average over the ellipsoidal profile to get the spherically 
symmetric one. On the data side, we hope to get the 
data to fill the gap in between 0 and 3 in future. As we 
see the model does not work very well, in practice it will 
be useful to come up with a parametrized form for the 
evolution of the profile. Also the paramertized form of 
halo can be used to improve the halo model. In the stan¬ 
dard halo model one assumes that all the matter exists 
within halos, and the virialized halo profile is used, such 
as the NFW profile for halos [30]. However at higher 


redshift, virialized halos are rare, and this assumption is 
not justified. One can improve the halo model using the 
proto-halo profile instead. 


III. BIAS WITH PROFILE CORRECTIONS 

We shall apply the fluid approximation to model the 
evolution of the dark matter and the galaxy field. The 
fluid approximation enable one to derive the nonlocal bias 
parameters [HEDIU] which results in better modelling 
of the halo power spectrum and bispectrum DJCnHET] 
and halo 3-point function [55] . In this paper we use halo 
and galaxy interchangeably. For dark matter, we will 
use the standard perturbation theory (SPT) results (see 
[5fi] for a review). In this framework, the evolution of 
the density contrast of the galaxy, i5g, and its velocity 
divergence 9^ are governed by the continuity equation 
and the Euler equation 


dr 


-b( 






d3fcid3A:2^D(k-ki2)a(ki,k2)0g(ki)^g(k2), 
d^kid^k2SB{k - ki2)/3(ki, k2)0g(ki),5g(k2), 


(13) 

(14) 




























where r is the conformal time, % is the conformal Hubble parameter din a/dr, ki 2 denotes ki + k 2 , and a and /? are 
the coupling kernels 


a(ki,k2) 


ki2 • ki 


k 


2 ’ 
1 


/3(ki,k2) 


■ k2 


(15) 


Here is the density parameter of matter. 

Eq. and are similar to the fluid equations widely adopted for modeling the evolution of multiple components 
m 12311571 158). except with the window function W, which is central to the results in this paper. We also note that 
in [15j . a similar modification of the Euler equation was proposed, in which the authors argued the forced for halos 
should be biased. However, the physical origin of this modification and its form are quite different from that in m- 
Another important difference from [TFl [2S] is that we do not impose the peak constraint in the evolution equations. 
The proto-halos after initial identification, they simply evolve following Eq. and 

The introduction of W is to model the fact that although 6g denotes the density contrast of the spatial distribution 
of the CM of the halos, each individual halo consists of a collection of particles. Thus the force on the CM of the halo 
should be the average force acting on all the individual particles in the halo. Hence in real space the effective source 
of the gravitational force for a finite-sized object is IT * instead of only S at the CM position of the object. In 
Fourier space, it is given by the product between W and S thanks to the convolution theorem. This window function 
describes the profile of the object. We will use the window function/profile studied in Sec. |h[ 

To simplify Eq. |13| and |14| further, we introduce the new time variable y = InD where D is the linear growth factor 
for the dark matter satisfying the equation 


(PD 


+ n^-ln^nmD = o. 

dr 2 


(16) 


We note that « flm, with / = d\n D/dlna is a very good approximation for the epoch that we are interested in 
[59]. Using this approximation Eq. 13 and |l4| can be written as 


dy ® . 

^ d^k^d^k^Sjoi^ - ki2)a(ki, k2)0;(ki)5g(k2), 

(17) 

d9g 1~ 3^^^, 

^d3fcid3fc2<5D(k-ki2)/3(ki,k2)0g(ki)4(k2), 

(18) 


where 9g denotes 9g/{—f'H). In the rest of the paper, we shall abuse the notation and simply use 9 and 9g to denote 
9 /{—/%) and 9g/{—f'H) respectively. 

In the following subsections, we will solve Eq. and 18 to linear and second order respectively to reveal the 
effects of the window function on the bias parameters. In |14| . the continuity and Euler equation of the galaxy field 
together with the other two equations for dark matter were written in a concise form, and hence a general perturbative 
solution was obtained using the transient formalism [39] , thanks to the fact that the coefficients of the equations are 
not explicitly time-dependent. However, W is time-dependent as we see in Sec. [TT] Here we will solve Eq. [TT] and [TS] 
directly. 


A. Linear biases 


and (5^^^ we get 


We start from the linearized version of Eq. and [TS] 



(19) 

dyOP + \oP = 

(20) 


where the superscript (1) emphasizes that the field is 
linear. We have suppressed the explicit k-dependence. 


1. Velocity 


We will work on the Euler equation first because it 
depends only on 9g^K Integrating Eq. 


20 


in terms of W 


( 21 ) 

As we mentioned before, in this paper, we use star 

*( 1 ) 

to denote a quantity at some initial time. Thus 9g^ ' is 
the velocity divergence of the galaxy at the initial time. 
For convenience, we define the function /„ as 

Iniy)= r dy'W{y')e^y\ (22) 

^0 

thus we have 

9^i\y)=9<^^e-^^+ I., (23) 
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z = 5. 1, 0 



FIG. 7. Evolution of the linear velocity bias with static 
window (dashed) and that with window given by SC (solid). 
Three redshifts are shown z = 5 (blue), 1 (red) and 0 (green). 
The case with static window is almost constant (thus some of 
the curves are covered by the blue curves), while the evolving 
window case gives decaying, but not negligible 6v. 


where J„ denotes the integral 

J„(y)= (29) 

do 

The advantage of introducing is that it gives at least 
order correction, thus it represents the genuine halo 
profile correction. The first two terms in Eq. are the 
velocity bias evolution obtained in m, and the last term 
is new, which arises from the halo profile. In the limit of 
large y, the profile correction does not vanish, instead 
tends to 1 + |e“T J3 (y). If we assume that the window 
is static, we get W. 

In Fig. we plot &v at z = 5, 1 and 0. To set the 
initial condition, b* at z* = 49, we borrow the results 
from peak theory [25] 

b:= (30) 

where Wq denotes the Gaussian window function and s„ 
is the spectral moment defined as 

s„ = dTT J P{k)Wl{kRG). (31) 


where we have used the SPT result 

5 ( 1 ) (y) =5^6*'. (24) 


Note that we have normalized the linear growth factor to 
be 1 at the initial time so that y* = 0. 

The linear velocity bias b^ then is given byQ 




= 6 : 


where we have used 


0(l)(y) 

3 

+2 

3(1) _ 


2 

(25) 


(26) 


Correspondingly the initial linear velocity bias 5* is de¬ 
fined as 


b 


* 

V 


q*(l) 


q(l) 


(27) 


In Fourier space, the window function approaches 1 at 
low /c, thus it is convenient to express the integral in 


Eq. 25 in terms of IT — 1. Hence we have instead 


by — 1 + {b* — l)e ^ + 2 ® 


■^My). (28) 


^ The linear bias parameters are defined differently from that in 
peak theory ESI, where a smoothing window function is divided 
by. For example, the defined here is equal b^Ws, where Ws is 
a smoothing window and only by is the called the velocity bias 
in US). Our more “direct” definition is closer to the standard 
treatment, where the window is not explicitly written down but 
its effects will be included in by. 


Let’s clarify the reason that we set the initial conditions 
using the peak theory even though we dispute about its 
prediction at late time. Various studies, e.g. [131145] show 
that the density and velocity cross power spectrum be¬ 
tween halo and matter in the Lagrangian space can be 
well fitted by the functional form motivated by the peak 
theory. Thus one may think that here we only use an es¬ 
tablished fact from simulations. However we argue that 
the subsequent evolution can be modelled by the simple 
fluid approximation augmented with the window func¬ 
tion. 

In the plot, as an example we will consider halo of mass 
2 X 10^^ We map the top-hat window size to the 

Gaussian window size using the relation Rq = Rth / V5- 
We compare the case when the window is static, with the 
window size given by the Lagrangian size and the case 
in which the window size is evolved by the SG model. 
The difference in treatment is only in J 3 in Eq. 


28 


We 

note that when the window is static, the resultant b^ is 
almost constant over time. However, when the window 
is evolved by the SC model, b^ decays over time, but it 
is not negligible at late time. 

Although we have used Eq. 30 as in the initial con¬ 
dition, the contributions from the initial scale-dependent 
part is small compared to the ones due to Jn, even if 
we have used 5* = 1, we find that the results are quite 
similar to that in Fig. This highlights that the scale 
dependence of 6v is mainly driven by the late time halo 
profile. 

In the literature, the velocity bias is often associated 
with k^ correction to both the density and velocity bi¬ 
ases, e.g. these corrections can be derived from the 
peak model [25|. In [23], Zel’dovich approximation was 
used to displace the peaks to the Eulerian space, and 
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M=6.9e + 13 


M = 2.0e+14 


M=6.1e + 14 


M=3.7e+15 



FIG. 8. The integrand |j3 | as a function of y for various values of k, obtained from simulation data (filled circles). The results 
from the SC (solid) and the static window (dashed) are also shown. The simulation data is from Oriana. 


they found that the velocity bias remains constant over 
time. The numerical measurement seemed to be in favour 
of the peak model result m- Here we show that taking 
into account that halos are composite objects there is sig¬ 
nificant A:^-correction to the velocity bias and it does not 
decay away over time. When the static window is used, 
we also find that reduces to W in the long term limit. 
However, when the evolving SC model is applied, the ve¬ 
locity bias is not constant, as can be seen from Fig. 
Even in the static window limit, our result {W) is still 
different from [24], which gets Wh'^ instead. The rea¬ 
son for this difference is that we do not impose the peak 
constraint in proto-halo evolution |25j . These differences 
can be used to differentiate these two models. 

Alternatively we can express the linear bias parameters 
in terms of the time derivative of the profile. Integrating 
the integral in Eq. by parts, bv can be written as 

bAy) = W{y) + [bt-Wme-iy 

- f dy'W (32) 
Jo 

where W denotes dW/dy. This form shows that there 
are two contributions to the correction from the pro¬ 
file, one from W and another from W. The contribution 
from W is simply the trivial smoothing. However, nu¬ 
merically performing derivatives on sparse data can lead 
to noisy results. Thus we will only use the form in terms 


of Jn, such as Eq. 

As the velocity bias is mainly generated by Js in 
Eq. to gain insight into which part of the integral 
of Ja contributes most, we plot the integrand of Js, js 

i|(y) = et^(w(y)-l) (33) 

in Fig. 1^ Again, in the range 0 < j/ < 3, there are no 
data available. We also show the prediction from the SC 
model. SC predicts that the contribution to the results in 
that range is small, while the contribution around y S 
is the largest. However, we note that the SC results of¬ 
ten overshoots in this range. We also show the results 
obtained with the static window. Static window approx¬ 
imation is good for y <2, but it overestimates the results 
for y ^ 3. 

As both the window function W and the time integral 
Ja contribute to leading k^ correction, we would like to 
compare the magnitude of these terms. In Fig. pj we 
compare contributions from |IF — 1| and 
using both the numerical results and the SC model. The 
results are for z = 0. The SC model gives quite good 
description of the results from the data. In particular, 
the value of |e“T|j 3 | from numerical window and the 
SC model agrees quite well. 

For the SC result for |IF — 1|, we have compared a 
few prescriptions for the size of the window. Although 
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M = 6.9e+13 M=2.0e+14 M=6.1e+14 M=3.7e+15 




FIG. 9. The comparison of the contribution to from \W — 1| (bine) and 
to correction in the velocity bias in Eq. 


2 IJal (red), which is the main contribution 
The resnlts from simulations (symbol) and SC (solid line) are shown. The 
simnlation data is from Oriana. For W — 1 from SC, we have used computed with Av = 500. The quantities are evaluated 
are z = 0. 


28 


at 2 = 0 the size has not yet collapsed exactly to zero, 
using such a value gives the magnitude oi W — 1 much 
smaller than the simulation results. We have tried using 
Eq. [^with Av = 200, 380 and 500, and Av = 500 gives 
the best agreement with simulations. In Fig. we have 
shown the results obtained using Av = 500. In passing, 
if we simply use Eq. which is strictly only for matter- 
dominated universe, we get the results very similar to 
those from Av = 380. We have cross-checked the results 
using Carmen, and they are consistent with those from 
Oriana. 

We note that recently there are reports of measure¬ 
ments of velocity bias at late time p^l7) . In [15], the 
momentum was measured, and they found that the model 
with 6pk arised due to the peak constraint seems to fit the 
data better at high redshift such as z = 20 than the evo¬ 
lution model without it. At such high redshifts, the effect 
of the profile evolution is small as can be seen from Fig.[^ 
Thus if confirmed, this would shows that the profile cor¬ 
rection would not be the dominant effect of the velocity 
bias seen in simulations. In the study of velocity 
bias was measured with sampling bias correction applied. 
They found that the velocity bias at fc ^ 0.08 Mpc~^ h is 
slightly positive, with bv ~ 1.01. In our model, velocity 
bias can only be negative in the mildly nonlinear regime. 
However, as the number density of halos is low and halos 
are more inhomogeneously distributed in Eulerian space. 


thus it is hard to get an accurate volume weighted mea¬ 
surement. It is not clear that these measurements are 
free of artifacts. Thus we will keep these in mind and 
hope to report our own numerical comparison in future. 


8. Density 


We now turn to the density bias. Plugging Eq. 21 into 


Eq. 19 we have 




(34) 


It is useful to note that 


r dy'e^y'h{y') = ^[e“^4(y) - h+biy)]- (35) 
Jo a 


Using Eq. 35 we can simplify Eq. 34 to 
(5W(y) =C(i) + 2C«(l-e-5)-3J«[e-«73-/i]. (36) 
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z=5. 1, 0 



FIG. 10. Evolution of the linear density bias with static win¬ 
dow (dashed) and evolving window given by the SC (solid). 
Three redshifts are shown: z = 5 (blue), 1 (red) and 0 (green). 


0=5, 1, 0 



FIG. 11. Same as Fig. |10[ except the scale independent 
initial conditions 6 * = bi, and b* = 1 are assumed. 


Thus the linear density bias is given by 


= b*ie-^ + 2b*e-^(l-e-^) 
-3e-^(e-^/| - Ji), 

where b^ is defined as 

j*(i) 

Ug- 


* — 




sr(l) 


(37) 


(38) 


Or in terms of J„ using = Jn + (e"^ — l)/n, we have 


bi{y) = l + {bl + 2 b: - 3)e-^ + 2(1 - 6 :)e-^ 


-b 3e"*'Ji - 3e"T j3. 


(39) 


The first line in Eq. 39 is the same as the time evolution 
of linear density bias obtained in [14], while the second 
line results from the halo profile correction. Unlike the 
decaying terms in the first line, they do not decay away. 


In the long term limit, bi reduces to l-|-3e ^ Ji—3e = J 3 . 
If we assume that the window is static, we get bi = W 
in the long term limit. 


In Fig. the evolution of the linear density bias is 
plotted. Again we use the form of the initial condition 
motivated by the peak theory [55] 


bl = {b, + b^k^)WG{kRG)- (40) 


We take Rg corresponding to halo of mass 2 x 
Instead of using the peak theory results, 
we take 6 ^ = 15.9 and = 40.0(Mpc/i“^)^, which 
are obatined from measurement of the initial cross power 
spectrum Both the results from the static window 
and SC evolving window are shown, however, the differ¬ 
ences are very small. 

The initial condition term {b\ -\- 26* — 3)e“*', especially 
due to 6 *, is important at low k. In fact, at low k, it 
gives the decay of linear bias [5T] . The term due solely 
to 6*, 2(1 — 6*)e“^^/^ is negligible in the whole range of 
k shown. The sum of the two profile correction terms, Ji 
and Ja gives small overall correction. That is also the 
reason why the static and evolving window gives almost 
identical results. The reason that the prohle correction 
term gives much more significant effect for 6 v than 
for 61 is that unlike the case of 6 v, 61 at late time is 
still dominated by the reminant effect of 6 * because the 
magnitude of 6 J is much larger than that of 6 *. 


Unlike the case b^, the magnitude of the scale- 
dependent part of 6 * in Eq. 40 is significant compared 
to other contributions. To highlight its effect, we plot 
the results when the initial bias is scale-independent, 
i.e. b: = 61 / and 6 * = 1 in Fig. 11 This plot shows 
that the bump in Fig. 


10 


_ around k ~ 0.7Mpc“ h is due 

to the large magnitude of the initial b/^k^ term. The low 
k plateau is due to b^, from initial condition and the scale- 
dependent transition comes from Ji and Ja terms. 
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B. Second order biases 


To second order, Eq. 17 and 18 become 


de, 


( 2 ) 


dy 


ar(2) p 

- J ki2)a(ki,k2)6»W(ki)(5W(k2), 

+ = yd'fcid"fc2^D(k-ki2)/?(ki,k2)0(iHki)0W(k2). 


( 2 ) ( 2 ) 

We shall solve Eq. 41 and 42 perturbatively to obtain 6g ' and 9g . 


(41) 

(42) 


1. Velocity 


(2) 

We will start from Eq. 42 to compute 9g first. Using the dark matter SPT result 


^(2)(j^) =g2y^(2)^ 

and Eq. we can integrate Eq. to get 

9 ^^\v) = 9 f-\-'^+ + J d^hd^k 2 Sj,(k--ki 2 )P(k,M) 

X e-« {2(1 - e-^)0*(i)(ki)0*(i)(k2) - ^[(e-5/|(k2) - /i(k2))0*(i)(ki)5l^)(k2) + (ki o ka)] 


+ \J ^2/'e ” 2 /3(j;',ki)/3(y',k2)(5^^^(ki)5i^^(k2)|. 


Replacing 9g^^'^ by b*si^'^ and extrapolating to the present time, we have 

(y) = +^<5(2)(y)e-^J|(y)+ y'd3fcid3MD(k~ki2)ife2(ki,k2)(5(i)(ki)^W(k2), 


where Kg 2 is given by 


ife2(ki,k2) = e-^/3(ki,k2){2(l - e- = )5:(ki)5:(k2) -^[{e-i- /i(k2))6:(ki) + (ki o k2)] 

+ l[ d,ye~'^Iiiy,ki)l2iy,k2)Y 


(43) 


(44) 


(45) 


(46) 


where (ki -f-)- k 2 ) is a shorthand for a similar term obtained with ki and k 2 interchanged. To be general, we allow 
the initial linear biases to be scale-dependent. Note that Kq 2 is already symmetric in ki and k 2 , so symmetrization 
is not required. 

q(2) 


In terms of J„, 9g’ can be expressed as 

9f\y) = 9<^^e-"^ + J d^hd^k26j,ik-k,2)}C92{ki,k2)S^^\'i^i)S''^\'i^2), 

where ICe2 is given by 

^fl2 = Tp + Tg + Tg, -\- Tjj + Tj. 

The five types of terms Tp, Tq, Tg*, Tjj, and Tj are given by 

TG = (l-e-^)G2(ki,k2), 

rF = ^e-'^JUk)F2(ki,k2), 

Z 2 

Tfc. = /3(ki,k2){2e-T(l - e-5)&;(ki)6:(k2) + e-^(e^ + 26"^ - 3)(6:(ki) + 6:(k2)) 
- 2(e-^ - 2e-^ + - ^e-^y [(Ja (k2) - Ji(k2))6:(ki) + (ki o k2)] }, 

Tjj = / 3 {ki,k 2 )e~^ f dye~^^Ji{y,ki)J^{y',k 2 ), 


(47) 

(48) 

(49) 

(50) 

(51) 

(52) 


Tj =/3(ki,k2)e ^ J dy'e “ 1)(>^|(y'.ki) + J 3 (y',k 2 )), 


(53) 
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where k = ki 2 , and F 2 and G 2 represent the coupling kernels 


^2(ki, k2) 


5 

7 


1 ,ki 


ki 



G2(ki,k2) 


3 

7 


1 

+ 2^ 


( 


^2 


+ 


ki 



(54) 


with ^ = ki • k 2 . 

As a cross-check, we pause to consider the limit b* = 1 and Jn = 0. Then ]Ce 2 reduces to (1 — e 
in this limit , thus Eq. 


47 


matter field. On the other hand, in the long term limit y 
from 6*, and we end up with G 2 -f Tp -|- Tjj + Tj. 

In Fig. 


)G 2 . Note that 
reduces to because the galaxy field reduces to the dark 
00 , the transient terms vanish, in particular those arising 


12 we show the evolution of the kernel lCg 2 - In this plot, we have used the same parameters as those in 


the previous section and have set ki = k 2 and fi = —1/2, which corresponds to the equilateral triangle configuration. 
We have compared the case with the static window and the one with SC evolving window, and find that the high k 
corrections are quite different. In particular, the magnitude of the static one decreases while the evolving one increases 
over time. When b* is assumed to be scale-independent instead, the results are similar to Fig. 12 thus we do not 
show it here. 

We now look at the individual components of ICg 2 in details in this example. At low k, the only non-vanishing 
component is Tq and it is (almost) constant for the reshifts shown. The term Tp gives negative k^ correction and 
its magnitude is large among all the high k correction terms. The term Tj and Tjj are of opposite signs, but the 
magnitude of Tj is slightly larger. In particular, as the leading correction from Tjj is of compared to Tj, it is 
unimportant for k < 0.6Mpc~^ h. Tp and Tj are the largest scale-dependent correction terms, they are of similar 
magnitude but of opposite signs. The term with 6*, Tg- gives small negative contribution, which is negligible compared 
to the other correction terms. In fact for k < 0.6Mpc“^ h, the kernel ICg 2 is well captured by the sum Tp -\-Tg + Tj. 
This is similar to 6v, for which the term solely due to b% is negligible at late time (even at z ~ 5), and the dominant 
correction term comes from the J„-term. 


2. Density 

q(i) /i(2) „„ j c(1) 


(55) 


We now compute Integrating Eq. 41 in terms of 9'^\ 9^\ and yields 

5 ( 2 ) = j*( 2 ) f f d^kid^k 26 r,{'k-'ki 2 )a{'ki,'k 2 ) [ dj/'6'W(ki)5W(k2). 

Jo J Jo 

Using Eq. 45 replacing dg^^^ by b*Si^'^ and extrapolating to the present time, we have 

- e-"-) - II - e-^yi2] + J d^kid^k25uik~ki2)P{'kiM) 

X e-2^{2(l - 2e-5 + e-^)5:(ki)6:(k2) + ^ (k2) - 2e-5 Ji(k2) -b Ii (k2))6:(ki) + (ki o k2) 


dy'ie 


'2 — P 2 


^J3(j/',ki)/|(2/',k2)}(5W(ki)5('Hk2). 


(56) 


Making use of Eq. and we can compute the second integral in Eq. to get 

J d^kid^k2SY){k - ki2)a(ki, k2) dy' 9 ^p {y', ki)( 5 (^)(y', ka) 

= I d\id^k2Su{k - ki2)e-2*'a(ki,k2){26:(ki)6t(k2)(l - ) + 2&:(ki)&:(k2)(l - 2e-« + e"^) 

-36t(k2)(e-5j3(ki)-/i(ki)) +36:(ki)(e-^/3(k2)-2e-^/i(k2)-f/i(k2)) +36:(k2)[(e-^-2e-^)J3(ki) 

+ 2/i(ki)-Ji(ki)] -^"dy'e-^J3(y',ki)[e-^/3(j/',k2)-/i(2/',k2)]}5W(ki)<lW(k2). (57) 

Therefore, we have 

< 5 ( 2 ) =< 5 *( 2 )+ 20 *( 2 )(i-e-|)- 3 j( 2 )[g-^j, J d^kid^k2Suik-ki2)Ks2{ki,k2)S^^Hki)6^^Hk2), (58) 
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where Ks2 is given by 

Ks 2 {ki,'k 2 ) = e-^y{ 2 ie-y - 2e-* + l)(/ 3 (ki,k2) + a(ki,k2))6:(ki)6:(k2) 

+ 2(1 - e-^) 5 :(ki) 6 t(k 2 )a(ki,k 2 ) + 3 [e-^/ 3 (k 2 ) - 2e-^/i(k2) +/i (k2)]6:(ki)(a(ki,k2) + ^(ki,k2)) 
+ 3 [(e-*' - 2 e-^)J 3 (ki) +2/i(ki) - /i (ki)] 5 :(k 2 )a(ki,k 2 ) - 3 [e-*/ 3 (ki) - /i(ki)]6t(k2)a(ki,k2) 

dy'e~^l 3 {y',ki) (e "5 - e"^)J 3 (y', k2)/?(ki, k2) + (^e "^/3 (y', k2) -k2)) a(ki, k2) 

(59) 

In terms of J„, can be expressed as 

42) =^*(2)+20*(2)(i_e-|) + J d^kid^k 2 SBik-ki 2 )ICs 2 (.ki,k 2 )S^^'>iki) 6 ^^\k 2 ), ( 60 ) 

where Ks2 is given by 

JCs2{ki,k2) = - 5 ) - 3 (e-^ J|(k) - J2(k))]F2(ki,k2) 

+ 2(e-^ - 2e-* + l)(a(ki,k2) +/3(ki,k2))5:(ki)6:(k2) + 2(1 - e-5)&;(ki)6t(k2)a(ki,k2) 

+ [ 2 (e^ - 3 + 3 e-^ - e"^) + 3 ( 6 -*^ J3 (ka) - 26 "^ Ji(k2) + Ji (k2))] 6:(ki)(a(ki, ka) + ; 9 (ki, k2)) 

+ [2(e*'-2e» + 2 e -5 -e-^) + 3 ((e-^- 2 e-*)j 3 (ki) + 2Ji(ki)-Ji(ki))]6;(k2)a(ki,k2) 


+ [e^ — 3 + 2e ^ — 3(e ^ J 3 (ki) — Ji(ki))] 6 )|(k 2 )a(ki, k 2 ) + Ai + A 2 + ^3 j”: 


(61) 


where A 2 , and Ag represent 


Ai = -^l dy'e " Ja(y',ki)[(e ^-e 2 ) (y'^ k2)/3(ki, k 2 ) + (e ” (y'^ k 2 ) - Ji(y', k 2 ))a(ki, k 2 )], (62) 

^ ./n ^ ^ ^ 


A2 = - 3 j dy'e ” | 2 (e 2-e 2 )(e ? - 1 ) J3 (y', ki)/ 3 (ki, k2) 

- e~^)Ji{y',ki) + {ey - e~^)Ji{y',k2) - {e^ - l)Ji{y',k2) a(ki,k2)|, 


-.3 ey 

^2 ~ ^ 
1 


^3 = ^e-«(e5 - 1)4[5(2 + )2a(ki, k 2 ) + 2(6*^ + + 10)/3(ki, ka)]. 


(63) 

(64) 


When the limit &* = 1, 6 * = 1 and J„ — 0 are taken, Sg^^ in Eq. 60 reduces to as the galaxy field becomes the 
dark matter field. In the large y limit, we have 


/CS2 = 


3 

110 


3e-2y(^e-^ J|(k) - J 2 (k))]F 2 (ki,k 2 ) + e-^y(Ai + A 2 + A 3 ). 


(65) 


As only the symmetric part of the kernel /Cs 2 contributes to the integral in Eq. |60[ we need to symmetrize the 
kernel as 


/C^2(ki, k2) = ^ (/C52(kl, k2) + /Cs2(k2, ki)) . 


( 66 ) 


However, to reduce the length of the formulas, we do not explicitly symmetrize them. However, in the final results, 
we always use the symmetrized kernel. 

Similar to that in Ref. |14j . we define the terms that deviate from the local biasing prescription as nonlocal terms. 
Thus at second order, the nonlocal terms are defined as 


= e> - ((>.«<"> + '!(«<■>)’). 


(67) 


At second order, the nonlocal terms are induced by the initial linear bias, and the initial second order biases do not 
generate new terms m- Eor convenience we consider 


5(2)_6i< 5(2) = j*(2)+20*(2)(i_e-i)+ J d^kid^k2Suik-k,2)x(kuk2)Siki)dik2) 


( 68 ) 
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where XS 2 is given by 


x(ki,k2) = e-^y[[ - - (bUk) + 26:(k) - 3)e^ - 2(l - bt{k))e^^ 

+ ^e- « - 5 - 3e- 5 J5 (k) + 3 J 2 (k) - Ji (k) + 3e« Ji (k)! Fa (ki, k 2 ) 

+ 2(e-^ - 26-5 + l)(a(ki,ka) + /3(ki,k2))6:(ki)6:(k2) + 2(1 - e-«)5:(ki)6t(k2)a(ki,k2) 

+ [2(e5 - 3 + 3e-^ - e"^) + 3(e-V3 (ka) - 2e-^ Ji(k2) + Ji (ka))] 6:(ki)(a(ki, ka) + /3(ki, ka)) 

+ [ 2 (e*'- 2 e* + 2 e-« -e-^) +3((e-*'-2e-5)J|(ki) + 2 Ji(ki)-Ji(ki))]&;(k 2 )a(ki,k 2 ) 

+ [^y — 3 + 2e 2 _ 3 ^e (k 2 )a(ki,ka) + Ai + A 2 + ^ 3 j"j (69) 


In Eq. 68 we have used bi given by Eq. 


39 


The initial second order biases are hidden in S, 


.*( 2 ) 


Suppose that the initial conditions are given by 


S*(A = ^-l(Si^)f + blsP, b: = l, 


(70) 


where b^ and b* are scale-independent. In other words, we suppose that the initial density biases are local in Lagrangian 
space, and there is no initial velocity bias. If we also neglect all the J„ terms, then Eq. [^is simplified substantially 
and we end up with m 




(71) 


where 62 and 72 are given by 


b2=ble-^y, 72 = ?(6t-l)e-2^^(e^-l), 


(72) 


and Q 2 denotes 


g2{k) = J d^hd^k2Su{k-k^2){f^^-l)S^^'>{k,)6^^\k2). 


(73) 


In particular, because there is no velocity bias as the dipole term vanishes. 

We now consider the correction to these results due to initial scale-dependent biases and the corrections arising 


from the profile corrections. For the initial conditions, we will assume that b* and b* are given by Eq. 30 and 40 
respectively, and for simplicity 63 is a constant and 

We first define the scale-dependent parameters of the initial biases 


el{k)=bUk)-b,, e:(fc) = &:(fc) - 1 . 


(74) 


Then we have 


42) _ 6,5(2) 


- I 2 Q 2 = y ki2)^/>(ki,k2)5(^i(ki)5i^i(k2), 


(75) 


where the kernel '0 is defined as 


^/i(ki, ka) = Te* -I- Tf. -|- Tj 


(76) 
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with various terms given by 


T,. = ie-22^[-e*'et(k)i^2(ki,k2) + (e"-l)et(k2)a(ki,k2)] +(ki ^k2), (77) 

T,. = ^e-^y[2{e-y -2e-i+ l)(a(ki, k2) + ^(ki, k2))(e:(ki) + e:(k2) + e:(ki)e:(k2)) 

+ 2(e^ - e*')e:(k)F2(ki, k2) + 2(1 - e-5)e*(ki)(6, + (k2))a(ki,k2) 

+ 2(e« -3 + 3e-« -e-^)e:(ki)(a(ki,k2)+/3(ki,k2)) 

+ 2(6*^ - 2e5 + 2e-5 - e-^)e:(k2)a(ki,k2)} + (ki o k2), (78) 

Tj = [ - 3e-5 J5 (k) + 3J2(k) - 3eVi(k) + 3e« Ja (k)] i^ 2 (ki, k 2 ) 

+ 3[e-^j3(k2) - 2e-* Ji(k2) + Ji (ka)] 6:(ki) (a(ki, k2) +/3(ki, ka)) 

+ 3[(e-^ - 2e-«)J|(ki) +2Ji(ki) - (ki)]6:(k2)a(ki,k2) 

- 3(e-^ J3(ki) - Ji(ki))6*(k2)a(ki,k2) + + ^ 2 } + (ki o k2). (79) 


In Fig. we show the evolution of the kernel tjj and its components at 2 ; = 5, 1 and 0 respectively. The parameters 
used are the same as those in the previous section and we have set ki = k 2 and /i = —1/2. Similar to the case of 61, 
there is no noticeable difference between the case with static window and the evolving one. All these contributions 
peaks around k ~ 0.7 — 0.9Mpc“^ h. In this case, Tg* and T). are of similar magnitude but opposite signs, so they 
roughly cancel each other. Thus the net contribution is mainly given by Tj. The overall contribution of ip decays 
over time. 

Again to highlight the effects of the initial scale-dependent bias, we show the case when the initial condition is 
scale-independent, i.e. = bi, and b* = 1 are assumed in Fig. 


14 


Fig. 


13 


We find that the bump around k 

is no longer present, instead there is smooth transition from k ~ 0.2 to IMpc”^ h. 


1 Mpc ^ h 


The second order kernel will contribtute to the tree level bispectrum. As a quick check of the importance of the 
correction term, we compare the kernel ip with the nonlocal term kernel 72 (m^ ~ 1) ™ Fig. 15 The parameters used 


are the same as the previous ones. From Fig. 15 we can see the k^ correction at low k. However, the k^ correction 


term starts to surpass the nonlocal term at fc ~ 0.4 Mpc ^ h. Hence one should also find the signature of the 
correction term in the halo bispectrum at fc > 0.2Mpc~^ h. 

I 


z=5. 1. 0 z=5, 1, 0 




FIG. 12. The evolution of the kernel of for the case with FIG. 13. The evolution of the kernel of ip for the case with 

static window (dashed) and the SC evolving one (solid). The static window (dashed) and the SC evolving one (solid). How- 

parameters fci = fca and /r = -1/2 are used. Three redshifts ever, these two cases are indistinguishable. The parameters 

are shown, 2 = 5 (blue), 1 (red) and 0 (green). are set such that fei = k 2 and ^ = -1/2. Three redshifts are 

shown, 2 = 5 (blue), 1 (red) and 0 (green). 
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z=5, 1, 0 



FIG. 14. Same as Fig. |13[ except with scale-independent 
initial conditions 6* = and 6* = 1. 


z = 5. 1, 0 



FIG. 15. The ratio between the kernel ip ^nd 72 “ l)i Ihe 
kernel of the nonlocal term C/ 2 . The parameters ki = k 2 and 
= —1/2 are used. Three redshifts are shown, z = 5 (blue), 
1 (red) and 0 (green). Solid line for positive value and dotted 
line for negative one. 


IV. CONCLUSIONS 

Recent measurements of the velocity bias suggest that 
the velocity bias of the halos is non-negligible at the 
weakly nonlinear regime k ~ 0.1Mpc~^/i at late time. 
In the time evolution model of the halo held, previously 
it was shown that it leads to decay of the initial veloc¬ 
ity bias using the point particle approximation so that it 
becomes negligible at late time. On the other hand, the 
peak model gives constant velocity bias over time. Thus 
the measurement seems to be in favour of the peak model 
result. 

It is often assumed that halos are point particles and 
focus only on their center of mass. Here we argue that as 


halos consist of a collection of particles, the force acting 
on its CM of the halo should be the force averaged over its 
constituent particles instead of only the force at the posi¬ 
tion of the CM. To take into account of the halo profile, 
we introduce a window function in the Euler equation. 
We find that the window function leads to non-negligible 

correction to the linear velocity bias. While the initial 

velocity bias decays away, the correction due to profile 
correction does not. In contrast, in the peak model, the 
imposition of the peak constraint leads to an extra con¬ 
stant scale-dependent bias. This difference can be used 
to distinguish these models. The profile correction also 
gives correction to the second order velocity kernel. 
For the density bias, the effect of the profile correction 
is not important at low k because the magnitude of the 
initial scale-dependent density bias is large in the peak 
model. Thus even at low z, the magnitude of the terms 
due to are the most important ones. Nonetheless, 
this implies that the k'^ correction is non-negligible for 
k > 0.2Mpc“^ h, especially for bispecturm. 

Since the window function is dynamical, we model it 
using the spherical collapse model. We also measure the 
evolution of the halo profile by constructing proto-halos 
at different redshifts. To our knowledge this is the first 
systematic numerical study of the evolution of the proto¬ 
halo profile. We hnd that the proto-halo profile evolves 
from a top-hat-like profile to an NEW profile. We find 
reasonable agreement between the spherical collapse and 
the numerical results. On the theory side, one may im¬ 
prove the modelling using ellipsoidal collapse model in¬ 
stead. Computationally, it would be useful to come up 
with a parameterization for the halo prohle at various 
epochs. 

Our work has highlighted the importance of halo pro¬ 
file and its evolution on bias. In theories such as the 
excursion set theory and peak model, window function 
are used to define halos in Lagrangian space. They are 
often assumed to be static and the window size given by 
the Lagrangian size even when they are transformed to 
the Eulerian space. The idea of profile evolution can be 
easily applied to these models as well. 

In our model, the effect of the window function cor¬ 
rection is most apparent in the velocity bias. Although 
there are some existing measurements of velocity bias, it 
is still hard because it is prone to sampling artifacts. We 
hope to report the comparison of our model with velocity 
bias measurement in future. 
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